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Abstract 

A kinetic flux-splitting procedure used in conjunction with locai ther- 
modynamic equilibrium in a finite volume allows us to investigate numer- 
ically discrete-velocity gas flows. The procedure, outlined for a general 
discrete-velocity gas, is used to simulate flows of the nine-velocity gas, a 
simple two dimensionai multiple-speed discrete-velocity gas, wherein a mul- 
tiplicity of speeds ensures nontrivial thermodynamics. After verifying the 
linear wave limit and the non-linear steepening of wavefronts, the stability 
and propagation of planar discontinuities in that model gas is studied. The 
supersonic-subsonic requirement for the stable propagation of a discontinu- 
ity, being kinematic in nature is the same in the model gas, as e.g., in a 
perfect gas. However, the finiteness of the velocity space in the model gas 
does not allow a translation of the above kinematic condition to the ther- 
modynamic requirement of increasing entropy across a compressive shock: 
a case of an entropy decreasing compressive shock in the model gas is pre- 
sented. Finally, the interaction of various types of waves — shock waves, 
rarefactions and contact surfaces — in the model gas are shown in a simula- 
tion of the shock tube problem. 

Introduction 

Regular discretization of the velocity space of particles constituting a 
gas results in a discrete-velocity gas. The computational evolution of such a 
system, especially the collision phase, is greatly simplified if the discretiza- 
tion is not elaborate. Lattice gases involve a discretization of the physical 
space as well, allowing them to be mapped directly on to a digitai com- 
puter, resulting in elegant computational models for fluids. The work by 
Frisch et. al.-*^, wherein the incompressible Navier-Stokes equations were 
recovered for the Frisch-Hasslacher-Pomeau (FHP) model, and the subse- 
quent success of using lattice gases'^'^ in simulating various fluid phenomena 
have renewed interest in these discrete-velocity models, first used in a flow 
situation by Broadwell'*'^. The predominant use of single-speed models to 
simulate incompressible flow notwithstanding, the macroscopic behavior of 
these models is compressible. Therefore, to investigate the applicability 
of discrete-velocity models to compressible flows, we have been studying a 
few multiple-speed models; the independent energy variable in the multiple- 
speed models allows nontrivial thermodynamics. We note that compressible 
flow dynamics of single speed models have been studied^" ^. 
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The computational technique with the attendant flux-splitting and 
total-variation-diminution^ ideas are outlined after briefly describing the 
nine-velocity gas-'^^"-'^^. Linear wave propagation and non-linear steepening 
and shallowing of planar waves — a familiar consequence of the convective 
nonlinearity in fluid equations — are then considered. Next, jump conditions 
associated with the model gas, time evolution of the jumps and how they 
relate to thermodynamics of a perfect gas are studied. We end with a sim- 
ulation of the shock-tube problem in the model gas, and a few conclusions. 

The Nine-Velocity Gas 

The nine-velocity gas, a two-dimensional model, consists of a large 
number of identical hard sphere (disk) particles, each of which takes on one 
of the nine allowed velocities 

Co = (0,0), Ci = (l,0), C2 = (l,l) 

C3 = (0,l), C4 = (-l,l), C5 = (-l,0) 

C6 = (-l,-l), C7 = (0,-l), C8 = (l,-1) 

Except for the simpliflcation of discretizing the velocity space, the model gas 
is much like a monatomic gas with a hard-sphere potential, and consists 
of free-flight interrupted by instantaneous collisions with other particles. 
The collisions, in addition to conserving mass, momentum, and energy 
individually, are closed under the velocity set. The four types of collisions 
possible are 

type 1: (cs, C7) ^(ci, C5), type 2: (04, cs) ^(c2, ce) 
type 3: (ci, C3) ^(cq, C2), type 4: (ci, ce) —^{05, cs) 
Collision type 3 is unique in that the precollision speeds are different from 
the postcollision speeds, and this provides the cruciai mechanism for equi- 
libration between the various particle speeds. 

If n [ni=n[ci), i = 0, 1, . . ., 8] represents the full velocity distribution 
function, then at thermodynamic equilibrium, from a detailed balancing of 
collisions, there are flve relations among the nine population densities: 

nius = non2, 713715 = noUi, n^wr = none 
nrni = nons, n^n^ = n^n-j 

From Eq. (1) only four of the ni are independent; a set of four such particle 
populations is denoted by m, where m is a Maxwellian state. 

The vector of mass, momentum, and energy, denoted by F, is given by 

F = ( X] X] X] "■"'^^ ) = 

\a=0 a a / 

2 

where E — e -\- ^ and G, the flux of F is given by 

G= ^?^aCa, ^riaCaCa, ^riaC^Ca (3) 

\ a a a / 
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With the above defìnitions of F and G, and the equilibrium relations 
Eq. (1), the Euler equations for the model gas may be written implicitly 
in terms of m as 

+ V • G(m) = (4) 
Pressure in the model gas, as described above is-'^'* 



4e 



correct to 0[u^): Note that at e = |, Eq. (5), the equation of state of 
the model gas reduces to that of an ideal gas. Because of the kinematic 
dependence of the thermodynamics, the speed of propagation of sound in 
such a gas is dependent on the macroscopic flow velocity, and because 
of that, there is not a unique speed of sound to which velocities can be 
referred-'^'*. This makes the definition of a Mach number in a discrete- 
velocity gas neither naturai nor unique. 



Simulation Scheme for a Discrete- Velocity Gas 



We briefly outline the computational technique here; the details can be 
found in Nadiga and Pullin®. For simplicity, a flow in one spatial dimension 
is considered: the flow is along the x axis, coincident with ci. Instead of 
dealing with individuai particles, the flowfield is tiled with a linear array of 
cells, with each celi interacting with its adjacent neighbors through a flux 
of F, the vector of mass, momentum, and energy. The cruciai step of the 
technique consists of splitting G, the flux of F, based on the signs of the 
discrete- velocities of the particles into G"'" and G~ and then computing the 
split fluxes in accordance with locai thermodynamic equilibrium. The time 
evolution then simply consists of noting that F is conserved in time and 
updating the particle distribution function in each of the cells accordingly. 
To illustrate the procedure, we write a first-order scheme based on the 
above ideas using an Euler time step: 

m(a;, t + Ai) = m(a;, t) 

[JFm]-' { G+ [n+{x, t)] - G- [n- {x + Ax, t)] (6) 
-G+[n+(a; - Ax,t)] + G-[n-(a;,i)] } 

JFm is the Jacobian of the transformation from F to m, and follows from 
the definition of F in Eq. (2), n"'" is the distribution function of particles with 
a positive u-velocity component and n~ that of particles with a negative 
u- velocity component. According to Eq. (6), the distribution function at x 
at time t + At is different from that at time i by 1) the departure of particles 
from X due to a non-zero u- velocity, terms 1 and 4; 2) the arrivai of particles 
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with a positive u-velocity from e — Ae, terni 3; and 3) the arrivai of particles 
with a negative u-velocity from x + Ae, terni 2. However, the important 
difference from the usuai discrete- velocity gas evolution is the fact that the 
arrivai and departure of the particles is in accordance with the equilibrium 
split-flux of F. In the present case of the model gas, with the flow along ci, 
the definition of G used in conjunction with Eq. (1) gives the expressions 
for the locai equilibrium fluxes, and while particles with velocities ci, C2, 
and cs contribute to 0+, particle types C4, C5, and ce contribute to G~. 
The minmod^ limiter is defined as the binary operator: 

minmodb, q) = sgn(p) | ° • n , , n !! ^^"J^ì ^ ^^"J^ì (7) 
° ^ ' [ min{|p|, if sgn(p) = sgn(g) 

where sgn(p) is the sign of p and |p| is the absolute value of p. In the present 
usage, p and q are the values of the slopes at the centroid of a celi: p being 
the backward slope and q the forward slope. The full minmod limiting 
procedure simply consists of applying the above binary operation to each of 
the cells at any given time step to obtain the slopes for any relevant quantity 
at each of the celi centroids. Using the minmod limiter to interpolate F 
gives rise to slope-limited schemes, whereas using it to interpolate G gives 
rise to flux-limited schemes. In either case, the minmod limiting strategy 
results in a diminution of the total variance of the interpolated field, giving a 
stable second-order scheme. There being no inherent advantage of either the 
flux-limited or the slope-limited scheme over one another, we use the flux- 
limited scheme. The Euler time stepping in Eq. (6) was purely illustrative; 
we use a fourth order Runge-Kutta scheme in ali the cases. 

Linear Wave Equation Limit 

The initial disturbance, imposed over an uniform state [po=l, eo=0.33, 
uo=0.00) of the model gas in a periodic [0,1] domain in the direction of 
Ci, is in the form of a hyperbolic secant variation in density. The time 
evolution of pressure, as defined in Eq. (5), is shown in Fig. 1. The behavior 
of the other quantities is identical, in concurrence with linearity of the 
(nondispersive) wave behavior. Throughout the computation, the values of 
mass, momentum, and energy integrated over the domain are ali conserved 
to better than one part in a million. On reversing time, the two waves 
recombine to give the initial disturbance correct to the irreversible diffusion 
arising out of the numerical viscosity"'^'^ of the technique, but which can be 
interpreted in terms of the viscosity of the model gas-'^'* (and imperfections 
in the initial data). 

Nonlinear Wave Steepening 

The best setup for this is of course a simple compression wave — only 
one of the three families of characteristics not parallel — propagating into 
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an uniform state. Not knowing the Riemann invariants in the model 
gas prevents such an explicit setup. Therefore, the time evolution of an 
initial hyperbolic tangent interface between two equilibria (which satisfy 
the Rankine-Hugoniot jump conditions-*^^ arising from Eq. (4)) is studied. 
The gradients of ali of the flow variables at the inlet and the outlet are 
prescribed to be zero. Figure 2 compares the steepness of the pressure 
profiles. The steepening of the front is clearly indicative of the nonlinearity: 
the speed of propagation of the disturbance varies as the magnitude of 
the disturbance (say p) [cf. d[u + a) = j^dp for a perfect gas]. Figure 2 also 
shows how an initially steep expansion front gets shallower with time. (The 
situation in neither of the above two cases is exactly that of a simple wave 
because ali that can be said is that the variation spans across two uniform 
regions.) 



Although ali of the discussion in this paper assumes the frame of refer- 
ence of the model, to write an expression for the jump across an infinitesi- 
mally thin discontinuity, occuring in Eq. (4) and moving with a velocity U 
in the direction of ci, we momentarily fix the frame of reference with the 
discontinuity; then d/dt can be replaced by —Ud/dx and integrating the 
resulting equations from ahead of the shock to behind it gives the shock 
jump conditions-*^^ 



[G(m)l = C^[F(m)l or equivalently [G(n)l = C^[F(n)l (8) 



The computational setup used in studying the discontinuities is ex- 
actly like the one used for studying the steepening of a compressive wave: 
the flux-limited, second-order scheme in the domain [0,1] with Neumann 
boundary conditions at inflow and outflow. The initial conditions in ali 
the cases is a step discontinuity Jq between the upstream state and the 
downstream state 



obtained as a solution of the algebraic system Eq. (8), unless stated oth- 
erwise. The subsequent time evolution of the discontinuity is studied. Ali 
of the quantities in this paper are in their nondimensional form: density 
corresponds to the average number of particles in a celi, velocities are non- 
dimensionalized by q, the unit speed of the model, and energies by . 



Jump Conditions 



where [e] = — x. 



Evolution of Planar Discontinuities 





Case 1 
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The initial condition in Fig. 3 

Jo : (1.000, 0.100, 0.000) — ^ (1.137, 0.152, 0.100) (10) 

has no explicit Information about the direction of propagation of the dis- 
continuity Jq or its speed. Jq satisfies the jump conditions of the Euler 
equations, Eq. (8), for ?7=0.828. Evolving the initial conditions, the jump 
Jo is seen to be preserved with time and it propagates with a speed ?7=0.828 
to the right, exactly as obtained from the model Euler equations. The prop- 
agating discontinuity is compressive: this follows from a calculation of the 
pressure given by Eq. (5) (the concurrence of the above definition of pres- 
sure with its mechanical counterpart has to be looked into) upstream and 
downstream. The entropy, defined using a Boltzmann-like H function as 

l 8 

s = y^niln(ni) (11) 



is seen to rise across the shock in this case. The characteristic veloci- 
ties upstream and downstream are Cu = (0.742, —0.742, 0.000) and Cd = 
(0.881, —0.554, 0.014), where they are expressed as (a;_|_, wq), the three 
characteristics of Eq. (4). 
Case 2 

Time evolution of the initial jump 

Jo : (1.000, 0.200, 0.000) — ^ (0.872, 0.141, -0.100) (12) 

shows that the jump is stable in the sense that it does not break up into 
different jumps, but it is seen to disintegrate with time: the slopes of 
ali the quantities decrease continuously with time. The jump conditions, 
Eq. (8), do not rule out such discontinuities; and, therefore, a possible ex- 
ternal constraint has to be invoked. Kinematically, this is seen to be the 
supersonic-subsonic condition: the discontinuity must travel faster than 
the speed of propagation of an infinitesimal disturbance in the direction 
of propagation of the discontinuity upstream and slower than the speed 
of propagation of a disturbance in the same direction downstream. Thus, 
for a C_|_ discontinuity, U should be such that (w_|_)^ ^ U < (a;_|_)^. We 
note parenthetically that whereas a similar relation would hold for a C_ 
shock (but with the absolute values of the speeds concerned), a Co fam- 
ily of shocks are also possible by replacing subscript + with subscript 0. 
A discontinuity thus propagates stably if the effect of the passage of the 
discontinuity at a point is to increase the speed of propagation of a small 
disturbance in the direction concerned there; for a stable C_|_ discontinuity, 
> 0. The disintegration of the jump Jo, Eq. (12), which satisfies the 
jump condition is then explicable as due to the violation of the supersonic- 
subsonic condition. The characteristic velocities upstream and downstream 
are C„ = (0.775, -0.775, 0.000) and Q = (0.531, -0.886, -0.013), and the 
initial shock velocity U = 0.680. 

or li^+Ì^O (13) 
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The stable shock in case 1 is seen to satisfy the supersonic-subsonic condi- 
tion. 

Case 3 

Jo : (1.000, 0.400, 0.000) — ^ (1.292, 0.455, 0.200) (14) 

satisfìes the jump conditions, Eq. (8), for ?7=0.885. Evolving the initial con- 
dition, Jo is seen to be preserved with time, see Fig. 4, propagating with a 
speed ?7=0.885 to the right, again in concurrence with the Euler equations. 
The main difference between the present case and case 1 is that the entropy 
of the downstream is lower than that of the upstream. If the definition of 
entropy used is correct, the behavior of entropy in the present case may ap- 
pear counter to what is expected from the second law of thermodynamics. 
To better understand this, we consider the behavior of entropy in the more 
familiar case of a perfect gas. The kinematic supersonic-subsonic require- 
ment, [u + a] > 0, (more generally valid for any normal-*-® fluid), can be 
related to the variation of thermodynamic quantities across the shock by 

r 

d(u + a) = — dp. (15) 
pa 

where F is the fundamental derivative of gasdynamics ^ ['^^j i 
acoustic impedance and a the speed of sound. Therefore, across a shock. 



«1= / - 



dp (16) 



With F = > for a perfect gas, the supersonic-subsonic requirement 
can be satisfied only with pd > Pu or through a compressive shock, and the 
following relation-*^® ensures a positive entropy jump across a compressive 
shock * 

Mfl = W + 0(n^) where n = J^ (17) 



However, for a perfect gas, the entropy jump can be easily expressed in 
terms of the upstream and downstream quantities 

Since the model fluid being studied is two dimensionai, 7, the ratio of 
specific heats, is 2. Therefore 



f Pu ed\ 
\pd Su ) 



.loclogj^^ (19) 
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The Rankine-Hugoniot relations across a compressive shock (M^ > 1) in a 
perfect gas are such that > 1, ensuring a positive entropy jump across 

it. °" 

Since the particles in the model gas are hard spheres, like in an ideal 
gas, we assume the model gas to be a normal fluid, i.e., F > 0. Although 
we do not have the counterpart of Eq. (15) in the model gas, relating d(a;_|_) 
to dp, we point out that if the relation is similar to Eq. (15), then F > 
implies the possibility of only compressive shocks, and that is what is com- 
putationally observed. We cannot use Eq. (17) in the context of the model 
gas because it uses conjugate variables, temperature of specific energy and 
pressure of density, and their definition is presently not clear due to the 
finiteness of the velocity space. 

In Fig. 5, the entropy jump across a series of shocks in the model gas 
is studied. The entropy jumps are calculated by two different methods: 
on the X axis is the entropy jump calculated using the statistical formu- 
lation, Eq. (11), and on the y axis is log, (^^), from Eq. (19). This 

is done for two series of shocks: In the first series (triangles), the up- 
stream is held Constant at [pu, e^, Uu) = (1.0, 0.4, 0.0) and the piston moved 
into it at velocities ranging from 0.01 to 0.5. In the second series (dia- 
monds), the piston is moved at a Constant velocity of 0.2 into an upstream 
{Pu,Su,Uu) = (1.0,0.3 < Cu < 0.6,0.0). The linear least-squares fit (the 
solid lines) in the two cases indicate the validity of Eq. (19) as an approx- 
imation for the model gas, or more generally the correspondence between 
the thermodynamics of the model gas to that of a perfect gas, in the regimes 
considered. Although the negativity of the entropy jump in case 3 is ascrib- 
able to the finiteness of the phase space of the model gas, it is necessary to 
work out the thermodynamics in detail to fully describe such behavior. 
Shock Tube Problem 

Finally, a simulation of the shock tube flow with the model gas is con- 
sidered. The initial condition corresponds to a high-density, high-pressure 
driver section and a low-density, low-pressure driven section, separated by 
a diaphragm which ruptures at i = 0. The initial jump is 

Jo : (1.000, 0.500, 0.000) — ^ (0.200, 0.200, 0.000) (20) 

The gas in the two regions is the same. The end walls are modeled by mirror 
image sites across their actual location, resulting in specular walls which 
do not permit flow through them and are adiabatic. The time evolution of 
density is shown using a grey-level coding on the x-t piane in Fig. 6. The 
qualitative nature of the interactions of the various types of waves — the 
shock wave, the rarefaction fan, and the contact surface — in that picture is 
the same as in a perfect gas. 

Conclusions 

Using a new scheme to simulate discrete-velocity gases, we were able 
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to study the one-dimensional gas dynamics of a simple multispeed discrete- 
velocity gas. The scheme, based on locai thermodynamic equilibrium, has a 
viscosity which is interpretable kinetically, due to which discontinuities are 
captured without any spurious oscillations. While the model gas behaves 
qualitatively like a perfect gas in some situations, it exhibits seemingly 
anamolous behavior in others. A better thermodynamic analysis of these 
finite phase space models is the key to understanding them. 
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Fig. 1 An initial pressure pulse splits into a C+ and a C_ disturbance 
(a), each of which propagates with little distortion (parallel character- 
istics). Reversing time at t=0.2, the C+ and the C_ waves recombine 
to reconstruct the originai initial condition (b), correct to the difFusion, 
verifying linear superposability of the disturbances. 



Fig. 2 A compressive wave steepens with time (a), but a steep ex- 
pansion wave becomes shallower with time (b). In both fìgures, the 
profìles have been offset by the distance they have traveled to make 
the comparison. 



Fig. 3 An unsteady compressive shock in the model gas. The entropy 
increases across the shock. 



Fig. 4 An entropy decreasing shock. 



Fig. 5 A comparison of the statistical defìnition of entropy in the model 
gas to the thermodynamic defìnition for a 2-D perfect gas used for the 
model gas. 



Fig. 6 A shock tube flow in the model gas. 
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